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Abstract 

In this paper we present the general problem of identifying if a nonlinear dynamic sys- 
tem has a chaotic behavior. If the answer is positive the system will be sensitive to small 
perturbations in the initial conditions which will imply that there is a chaotic attractor in 
its state space. A particular problem would be that of identifying a chaotic oscillator. We 
present an example of three well known different chaotic oscillators where we have knowledge 
of the equations that govern the dynamical systems and from there we can obtain the corre- 
sponding time series. In a similar example we assume that we only know the time series and, 
finally, in another example we have to take measurements in the Chua’s circuit to obtain 
sample points of the time series. With the knowledge about the time series the phase plane 
portraits are plotted and from them, by visual inspection, it is concluded whether or not 
the system is chaotic. This method has the problem of uncertainty and subjectivity and for 
that reason a different approach is needed. A quantitative approach is the computation of 
the Lyapunov exponents. We describe several methods for obtaining them and apply a little 
known method of artificial neural networks to the different examples mentioned above. We 
end the paper discussing the importance of the Lyapunov exponents in the interpretation of 
the dynamic behavior of biological neurons and biological neural networks. 


1 Introduction 

In the companion paper [1] we described some findings about biological oscillators that have been 
presented in the recent literature. We also showed some examples of oscillator models. Here we 
want to review some models of chaotic oscillators with the goal of extending the analysis to time 
series (trains of action potentials) coming from biological oscillators where there are some hints 
that they are chaotic. 

There are many experimental situations where there is no idea of what the mathematical 
model of a system could be or where the form of the equations is known but the parameters are 
unknown. There is an extensive literature about methods for systems identification but they are 
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usually limited to linear models. Since the conditions for chaotic behavior arise from the presence 
of nonlinear elements, the use of linear methods is limited. 

In recent years there have been many advances in the understanding of nonlinear dynamic 
systems and that has produced many methods for identifying whether or not a given system 
is chaotic. For testing these methods various simple chaotic systems have been discovered or 
invented. In some of them the equations that govern the system are well known but in others 
only some sort of approximation is known. In the former case, it is possible to generate the 
corresponding time series with very high approximation, in the latter case, the measurements 
yield a sampled version of the corresponding time series. Having at hand the equations and 
the time series, or at least the time series, it is possible, in very different ways, to compute the 
asymptotic properties of the system. Two measures are used for this: the Lyapunov exponents 
in which a positive one indicates chaotic dynamics, and the attractor’s topological dimension 
which indicates topological characteristics and is directly related to the number of non- negative 
Lyapunov exponents [2], 

It is usual to ascertainl the existence of chaotic dynamic by means of visual inspection of the 
phase plane portrait. However, such method presents a considerable amount of^ uncertainty and 
subjectiveness. Taking that into account, it is important to have a quantitative method like the 
one provided by the computation of the Lyapunov exponents. 


2 Lyapunov Exponents 

To determine if a system possess chaotic dynamics it is necessary to know if it is sensitive to 
small perturbations on the initial conditions. When this occurs it is then impossible to predict the 
final state of the system after a finite time. To be able of characterizing a chaotic attractor it is 
necessary to establish quantitative measures concerning the sensitivity to initial conditions. The 
spectrum of Lyapunov exponents gives a method of quantifying the dynamics. The Lyapunov 
exponents describe the average rate of growing or shrinking of small perturbations in different 
directions in the state space. When the attractor has at least one positive exponent then it has 
the property of being sensitive to the initial conditions and it is called a chaotic attractor. 

There are several methods for computing the Lyapunov exponents. Wolf et al. [3] were the 
first in suggesting a method to compute them directly from the time series, without knowing the 
equations that govern the system’s dynamics. Kurths and Herzel [4] proposed another algorithm. 
However, in these algorithms the estimations are sensitive to the number of observations, to the 
sampling frequency and to the noise in the observations [3]. Trying to avoid these problems, 
Gencay and Dechert [5] designed an algorithm that computes the m Lyapunov exponents from 
an unknown m-dimensional dynamic system directly from a few observations on the attractor, 
in such a way that the estimation is robust even for certain amount of noise. This algorithm is 
based on a result by Hornik et al. [6] in which they show that the m Lyapunov exponents of a 
diffeomorphism that is 1 topologically conjugate to the process that generates the data, are also 
the m Lyapunov exponents of that process. To obtain a robust estimation considering both few 
observations on the attractor and the presence of noise, Gencay and Dechert [5] applied artificial 
neural networks with a cascade architecture. Such procedure is a non-parametric estimation that 
Hornik et al. [6] [7] have shown to be universal approximators, that is, they can asymptotically 
approximate a function and its derivatives. 
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3 Computation of the Lyapunov Exponents 

The Lyapunov exponents are constants, except for a zero-measure set, and describe the direction 
of nearby paths that converge or diverge in the state space of a dynamic system. The Lyapunov 
exponents A, are defined as the logarithm of the eigenvalues m V i = 1 , 2, • • • , m of the symmetric 
positive matrix 

A x = lim [y(x;f)‘ r T(x;f)l 1/2 \ (1) 

t— ► OO L J 

where the matrix Y is dependent on the differential equation that characterizes the dynamical 
system. A direct application of the above definition is not practical since the Y matrix grows 
exponentially due to the fast convergence of the columns in the direction of greater expansion. 
Using topological properties and an appropriate QR decomposition, the Lyapunov exponents are 
found by computing 



where R lt are the diagonal elements of the triangular matrix R. 

An alternative to the aforementioned algorithm is the use of neural networks which are capable 
of recovering a nonlinear map from a time series of iterates [8]. Here an unknown function is 
estimated and then it is possible to compute the Lyapunov exponents using the properties of the 
dynamic system [5]. 




FIG. 1. Phase Plane Portraits for the Logistic Map, the Henon Map, and the 
Lorenz System. 
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4 Examples of chaotic oscillators 

To be able of testing the effectiveness of the Lyapunov exponents one has to have at hand dy- 
namical systems with proved chaotic behavior. Many mathematical model systems are known, 
for instance: Henon, Rossler-chaos, Lorenz, Rossler- hyperchaos, Mackey-Glass, and others [3]. 
Physical model systems ar|e more difficult to produce, however we have the Belousov- Zhabotinsky 
chemical reaction [9] and Chua’s nonlinear circuit family [10]. Obviously, there We real physical 
chaotic systems but they are extremely complex as it is the case of atmospheric turbulence. 

TABLE I. True Lyapunov Exponents, equations representing the chaotic systems, 
initial conditions and parameters. 


True Lyapunov Exponents 

Logistic map 

Henon map 

Lorenz system 

0.673 

0.440 

1.51 


-1.620 

0.00 

-22.5 
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a = 1.4, b — 0.3 

a = 16.0, b = 45.92, c = 4.0 


TABLE II. Estimated Lyapunov Exponents. Logistic Map (q = 5, T = 100). 
Henon Map ( q = 10, T = 200). Lorenz System ( q = 15, T = 1000). The error is less 
than lxlO -3 . The non-spurious Lyapunov exponents are shown in boldface. 


Estimated Lyapunov Exponents 

P 

Logistic map 

Henon map 

Lorenz system 

1 

0.6794 



2 

0.6401 

0.3670 



-6.7823 

-1.5673 


3 

0.6350 

0.4502 

1.7285 


-2.4378 

-1.7331 

0.0411 


-2.4930 

-2.8164 

-23.72 

4 

0.6434 

0.4119 

1.5910 


-1.61106 

-1.4803 

-0.0710 


-1.7342 

-3.3658 

-20.973 


-5.0200 

-5.2263 

-80.325 

5 

0.6790 

0.4385 

1.4799 


-0.9073 

-1.5473 

0.0067 


-1.3544 ' 

-1.4859 

-20.977 


-1.8468 

-1.7651 

-60.702 


-3.2313 

-2.5605 

-92.584 
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We divided this section of examples in three parts. In the first part the ordinary differential 
equations for the Lorenz model were integrated using the SDRIV2 subroutines form Kahaner et 
al. [12]. For the QR decomposition use used the subroutines in [13]. In the second and third parts 
we computed the Lyapunov exponents by means of the neural networks approach [5]. 

For the first part of the section we show the results of computing the Lyapunov exponents 
for mathematical model systems with well known chaotic dynamics [3] [5] [8]. We computed the 
Lyapunov exponents for the Logistic Map, for the Henon map and for the Lorenz system [3] [5]. 
The corresponding phase plane portraits are shown in FIG. 1, and the equations, parameters and 
computed true values of the Lyapunov exponents are shown in TABLE I above. 

For the second part of the section we show the results when we assume that for the same 
chaotic systems than above, the equations are not known, only the time series. In TABLE II we 
show the computed Lyapunov exponents where the presence of one positive exponent indicates 
that the system is chaotic. 

According to the established notation for neural network architectures [11], p represents the 
number of nodes in the input layer and q represents the number of nodes in the hidden layer. 
The output layer has one node. The error is the quadratic average summation of the differences 
between the real and the estimated values for the time series, being T the total number of sample 
points in the sequence. 

For the third part of the section we show the results obtained when we used a nonlinear circuit 
of the Chua’s family with the parameters, components and initial conditions shown in FIG. 2. The 
temporal series was acquired by means of a digital storage scope, the x-coordinate is the voltage 
in the linear capacitor Cl and the y- coordinate is the voltage in the linear capacitor C2. A part 
of the phase plane portrait is also shown in FIG. 2, from it the temporal series was obtained using 
a sampling frequency of 500 Hz. 



CHUAS CIRCUIT 



FIG. 2. Nonlinear Circuit of the Chua’s Family and a part of its phase plane 
portrait that was plotted using as state variables the voltages in the capacitors Cl and 
C2. The phase plane portrait changes when the parameters in the circuit are varied 
between the limits denoted in the diagram. 

The estimated Lyapunov exponents for the Chua’s circuit are shown in TABLE III, when using 
the estimation for q = 15 and T = 2500. The error was less than 5xl0 -2 . Notice that Aj is a 
positive exponent which means that the dynamic behavior of the circuit is chaotic, as expected. 
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TABLE III. Estimated Lyapunov Exponents for the Chua’s Circuit. 


Estimated Lyapunov Exponents 

P 

Chua’s Circuit 

3 

Xi =+4.07216 
A 2 =-0. 32160 
A 3 =— 3.58462 


5 Concluding Remarks 

We have shown well known examples of mathematical models of chaotic attractors: Logistic, 
Henon, and Lorenz. We also showed an electronic model of a chaotic attractor: a circuit of 
the Chua’s family. In all these examples we computed the Lyapunov exponents as a measure of 
the system’s sensitivity to small perturbations in the initial conditions. For the example where 
we know the equations we used the standard method; for the other two examples, we used the 
method of Gencay and Dechert [5] that applies a neuronal network algorithm. We went from the 
easier examples to the difficult one, that is, the Chua’s circuit where the time series is obtained 
from direct measurements in the circuit. In the former examples the true values of the Lyapunov 
exponents are known whereas in the Chua’s circuit the Lyapunov exponents are estimated and 
to this has to be added the differences or variations, for any reason, in the parameters for the 
circuit’s components. The justification for such trouble is that in a real chaotic system there is a 
complexity even worse that in the electronic model. Therefore, the circuit provides a very valuable 
experience that afterwards can benefit the understanding of the real chaotic attractor. 

Recent experimental evidence points to biological neurons and biological neural networks as 
very likely sources of chaotic attractors. As in the case of chaotic chemical reactions the biological 
significance given to such behavior is speculative [9] [14]. Nevertheless, the application of the tech- 
niques described in this paper might be very useful for interpreting data from neurophysiological 
experiments where the electrical activity from many neurons is recorded simultaneously. 

Similar to Chua’s circuit situation, due to the usual difficult conditions of neurophysiological 
experimentation, in a biological neural network we can only obtain a short duration record of the 
compound time series (train of action potentials). The individual time series have to be separated 
and that procedure produces and additional source of uncertainty and every tool available for 
interpreting the results is welcomed [15]. 

As a very simple example (for details see companion paper [1]) let us consider the rhythmic 
firing single neuron (1) shown in the upper plot in FIG. 3. This time series was obtained from the 
simulation of a mathematical model for a single neuron and two of its phase plane portraits for 
two different physiological parameters (V — /, V — h) chosen as state variables show very clearly 
the possibility of chaotic behavior. In the lower plot in FIG. 3 we show another simulation of 
the rhythmic firing of a neuron (2) that belongs to a recurrent- ring network composed of single 
neurons like the ones given in ( 1 ). From the two phase plane portraits we cannot conclude that 
the neuron (2) in the network is chaotic and the conclusion about neuron (1) required of an 
expert. This ambiguity can be surmounted if the Lyapunov exponents are computed for these 
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time series and from there the importance of making such calculations and having experimental 
testing circuits for improving the confidence in the results. 



FIG. 3. (1) Rhytmically firing single neuron. (2) Rhytmically firing neuron be 
longing to a network formed with type (l) neurons. 

On the other hand, in biological experiments there are problems similar to the ones present 
when making measurements in the Chua’s circuit. When doing an extracellular recording, the 
duration of it is limited to a few minutes and afterwards a considerable amount of preprocessing is 
required to get to the individual contribution of each neuron recorded [16] [17]. The calculation of 
the Lyapunov exponents for these individual contributions could be added to help understanding 
the functional role of oscillatory neurons and oscillatory networks. That is the work that we are 
about to pursue. 
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